Resolving Cosmic Neutrino Structure: A Hybrid 
Neutrino TV-body Scheme 



Jacob Brandbyge 1 , Steen Hannestad 1 

E-mail: jacobb@phys.au.dk, sth@phys.au.dk 

department of Physics and Astronomy, University of Aarhus, Ny Munkegade, 
DK-8000 Aarhus C, Denmark 



Abstract. We present the first simulation capable of resolving the structure of 
neutrino clustering on Mpc scales. The method combines grid- and particle-based 
methods and achieves very good accuracy on both small and large scales, while keeping 
CPU consumption under control. Such simulations are not only ideal for calculating 
the non-linear matter power spectrum but also particularly relevant for studies of how 
neutrinos cluster in galaxy- or cluster-sized halos. We perform the largest neutrino N- 
body simulation to date, effectively containing 10 different neutrino hot dark matter 
components with different thermal properties. 
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1. Introduction 

Neutrinos are among the most abundant particles in our Universe, and since at 
least two of the three neutrino mass eigenstates have masses much larger than the 
current temperature, neutrinos contribute to the matter density and are important for 
cosmological structure formation. Oscillation experiments have established two mass 
differences, Am^ ~ 7.6 x 10~ 5 eV 2 and (Am^l ~ 2.4 x 10~ 3 eV 2 [TJ, from which the 
heaviest mass eigenstate must have a mass of at least 0.05 eV. Even such a relatively 
small mass will have the effect of suppressing the power spectrum of matter fluctuations 
by ~ 5%, an effect which is significantly larger than the precision with which the matter 
power spectrum can be measured in upcoming surveys. 

Consequently this has led to a much increased interest in gaining a detailed 
understanding of how neutrinos affect structure formation. In linear theory the effect 
is extremely well understood. However, these results apply only on very large scales, 
k <C 0.1 WVlpc -1 , whereas most of the cosmologically relevant information from large- 
scale structure surveys is in the range k ~ 0.1 — 0.5 /iMpc -1 . On these intermediate 
scales it is mandatory to correct for non-linear effects, even if surveys are carried out 
at intermediate redshifts where non-linearity is weaker. In the semi-linear regime the 
power spectrum can be estimated using analytic methods (see [21 EH II])- 

However, in the fully non-linear regime the most accurate, but also most time- 
consuming method is to use full iV-body simulations with neutrinos included in a self- 
consistent way. In a previous publication we carried out a large suite of simulations with 
neutrinos included [5] , and found a significant non-linear correction caused by neutrinos. 

In a recent paper [H] we showed that for realistic neutrino masses, quantities such 
as the power spectrum can be calculated very reliably on all scales, using a grid-based 
method which tracks the linear neutrino density contrast on a grid while using the full 
non- linear structure of CDM. 

The problem with both the particle-based and the grid-based methods is that they 
cannot reliably probe the small-scale structure of neutrinos. In particle-based codes 
the problem is that the thermal velocity of the neutrino component is so large that 
it introduces noise at an unacceptable level. The grid-based method does not resolve 
neutrino bound structures and therefore by construction does not allow us to probe for 
example the neutrino content of a halo. 

In the present paper we present a hybrid method which retains the good features of 
both methods and still runs at an acceptable speed. The idea is to start with neutrinos 
described on a grid, and then convert part of the grid to iV-body particles when the 
thermal motion of neutrinos decreases to a few times the flow velocities in the simulation. 
The neutrino iV-body particles are created in different momentum bins, with individual 
transfer functions and thermal velocities. We therefore present iV-body simulations 
with 15 different neutrino hot dark matter species. In a following paper we will use this 
hybrid method to investigate neutrino clustering on sub Mpc scales. 

In section [2] we outline the theoretical framework and problems with combining 
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particle and fluid approaches. Section [3] gives a technical description of the 
implementation of the hybrid method and Section H] presents our results. Finally, 
Section [5] contains our discussion and conclusions. 

2. Theory 

2.1. The Boltzmann equation 

In this subsection we will briefly outline how the evolution of massive neutrinos is 
followed in linear theory. The notation is identical to [7]. We use the metric in the 
conformal Newtonian gauge 

ds 2 = -a 2 (l + 2^)dr 2 + a 2 (l - 2<j))dx 2 . (1) 

In a perturbed universe the phase-space distribution function is expanded to first 
order as follows 

/ = /o + §^T = / (l + \EO, (2) 

with the perturbation parametrised by \l/ = — dlnfo/dlnq 5T/T, and the zeroth-order 
Fermi-Dirac phase-space distribution function given by 

^ = ^TT' (3) 

Here q\ = api, where pi is the proper redshifting momentum and T is the neutrino 
temperature today. 

After neutrino decoupling the collisionless Boltzmann equation evolves the 
distribution function as 

dj_ = df_ dx^_ df_ dq_ df_ dni_ df_ = Q 
dr dr dr dx l dr dq dr dni 

Expanding the perturbation $ in a Legendre series the perturbed neutrino energy 
density is given by a weighed sum over neutrino momentum states 

5p v {k) = 4vra- 4 / q 2 dqef ^ , (5) 



where e = (q 2 + c^m 2 ) 1 / 2 . Note that /o and not / is used to weigh the individual 
momentum bins, an assumption which breaks down when the gravitational flow velocity 
approaches the thermal velocity of the individual momentum states. 

From the Boltzmann equation the \!Vs are related to each other and the metric 
potentials by 

iff qk fiTf 2 uf ^ ek i dln f° 

— — - , I > 2. (8) 

1 e \2l- 1 2/ + 3 + V' - v ; 
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The second term on the right hand side of the equation for \&o encodes the effect 
of structure forming in evolving gravitational potentials, whereas the first term 
incorporates the effect of velocity on structure formation. The change in velocity is 
affected by the three terms on the right hand side of the equation for vE^. The first term 
encodes the effect of flow velocities redshifting in an expanding Universe, whereas the 
second term, found from the hierarchy of ^j's, incorporates the effect of momentum (and 
redshift) dependent neutrino free-streaming. These two terms give rise to less structure, 
whereas the last term in the equation for gives the acceleration as a gradient of the 
gravitational potential. 

The quantity calculated in linear theory is the transfer function (TF) which evolves 
the initial primordial perturbation ^ as 

Mk,Q,*)=T(k,q,z)* I (k,q), (9) 

with \1/q given by 

With adiabatic initial conditions 5T/T = <5„/4 and 5 c d m = 5b = \$v = f&y 0- 



2.2. Qualitative behaviour 

Eq. (jlOl) states that neutrinos with different initial momenta have different initial 
perturbations, also on large scales. Specifically, W Q — »■ for q/T — > 0, meaning that 
neutrinos with zero (read: almost zero) momentum are unperturbed. The reason for 
this is that at each spatial location the neutrino distribution is initially locally a Fermi- 
Dirac distribution with zero chemical potential, as is also reflected in \1/q. As q/T — > 
the Fermi-Dirac distribution always tends to / — > 1/2, independent of T, and this means 
that the distribution function is unperturbed. 

According to the Boltzmann equation, neutrinos with q = always stay 
unperturbed. When gravity acts on neutrino particles with q — 0, thereby perturbing 
them, their momenta change, and they are therefore shifted out of the q = momentum 
bin, leaving this bin unperturbed. This fact is contained in the Boltzmann hierarchy for 
the \£Vs through the term d\nfo/dlnq, which multiplies the gravitational potentials (<ft, 
ip) so that different neutrino momentum bins are affected differently by gravity 0- 

The effect of linear as well as non-linear forces on the momenta of neutrino particles 
changes the shape of fo beyond linear order. Therefore the weighing of the individual 
\EV S in the momentum integral changes. This means that the value of 5p changes, 
not only at non-linear wavenumbers but also in the linear regime. Therefore the linear 
theory average neutrino TF is not even valid at very large scales, though becoming more 
accurate for lower neutrino masses. 

| Note that this is purely a phase-space phenomenon and has nothing to do with the gravitational 
force felt by particles. 
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2.3. Converting the density grid to N-body particles 

When we convert densities on the grid into particles we are faced with a potential 
problem. Imagine that neutrino energy in the first % momentum bins are converted 
to particles at some redshift, while keeping the remaining neutrinos on the grid. At 
some later time, some of the created particles will have been moved to higher velocities 
by gravity and therefore really belong to the grid, whereas some density on the grid 
would have moved to lower velocity, and should have been converted to particles. This 
"leakage" is problematic because the evolution of the grid is done in CAMB [8] while 
the evolution of particles is followed in gadget-2 [5]. 

The problem is even more severe because the leakage to higher and lower momentum 
bins are qualitatively different: Imagine we dump all neutrino grid bins with momenta 
less than g cut to A-body particles at a certain redshift, and let the particles evolve. Due 
to energy conservation, particles gaining energy above g cut are falling into gravitational 
potential wells, and are typically clustered, while particles loosing energy are moving 
out of potential wells, and they are typically much more homogeneous as a population. 
Hence one cannot simply artificially decrease the velocities of the high energy particles 
to account for the leakage across g cut between grid and particles, because the topology, 
or specific position in phase-space, of particles crossing the momentum boundary from 
above and below is very different. 

Ideally, the neutrino A-body particles should only sample momentum space up 
to q cut . In the A-body simulation neutrinos with momenta larger than q cnt should 
therefore not contribute to the gravitational source term (to avoid double counting). 
This can surely be done, by eg. zeroing the mass of these neutrino A-body particles. 
Furthermore, we get a zero counting of the particles which, if they had been created 
from the grid at q > q CVLt , would have moved to q < q cut . One could include zero mass 
ghost particles with initial momentum states covered by the grid as tracer particles 
in the A-body volume, and then turn their mass on if they leak down to momentum 
states smaller than q cut . But in practise the linear assumption of a constant / ~ f is 
violated beyond linear order in the A-body simulation: There is a significant net flow 
of neutrino mass to higher momentum states. Therefore, turning the neutrino mass 
on / off as particles cross q cut will lead to a net reduction of overall neutrino mass. 
Furthermore, particles moving to q > q cut carry large non-linear corrections which are 
not found in the corresponding linear grid bins, making the idea with ghost particles 
less appealing. 

However, a simple solution which works well in practise is to convert a sufficiently 
large number of momentum bins at the same time because that keeps the leakage 
between converted and non-converted bins to a minimum. We will discuss this issue 
in the next sections and demonstrate that errors can be kept under control at the 
required precision. 
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Figure 1. Momentum dependent neutrino TFs at z = for ^2,m v = 1.2 eV. 
3. Implementation of the Hybrid Neutrino Method 

We have used the following parameters in a flat cosmology: (fl c + Q u , Q^, Q\, h, erg, n s ) = 
(0.25, 0.05, 0.7, 0.7, 0.878, 1). The value for a 8 is for a pure ACDM model without 
massive neutrinos. We use values of ^2 m v = 0.6 eV and 1.2 eV with particular emphasis 
on the latter value, since a high neutrino mass affect structure formation more, thereby 
testing the accuracy of the hybrid method better. 

The primordial density perturbations are followed with the linearised Einstein and 
Boltzmann equations, solved in CAMB [8]. At z = 49 the CDM component is taken 
out of CAMB and followed in gadget-2 [9] run in the TreePM mode. The A^-body 
initial conditions are generated with the Zel'dovich Approximation [TU] as well as a 
second-order correction term [XT] added to the CDM particles only. 

As input to the A^-body code the linear neutrino TF is needed. Normally this 
TF is found from Eq. ([5]) by integrating over all momenta. However, in our case we 
need the momentum dependent TFs because when particles are extracted from the 
grid they are taken in a certain momentum range and the neutrino TFs are highly 
momentum dependent (because the free-streaming length depends on q) . In Fig. [T] we 
show the momentum dependent TFs for Yl m f = 1-2 eV at z — 0, from which the 
neutrino momentum dependence on the power spectrum is clearly seen. The effect of 
free-streaming at small scales is most pronounced for the highest neutrino momentum 
states, but these states do have the largest TF at large scales. 

To follow the neutrino component we combine the neutrino particle and grid 
representations, first explored in [3] and [BJ, respectively. Initially, the neutrinos are 
represented on a grid, and there are no neutrino A^-body particles. When the neutrino 
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thermal velocity falls below a few times the average gravitational flow velocity, calculated 
from the CDM iV-body particles, the lowest momentum part of the neutrino grid energy 
is converted into iV-body particles. Here there are 3 things to consider. 

First, the energy on the neutrino grid should be decreased. We have 15 neutrino 
TF bins (the standard output from CAMB, the number can easily be modified), where 
q/T runs from 1 to 15 (upper limits) in integer steps, corresponding to a bin width of 
~ 125km/s in comoving coordinates for '^2 i m y = 1.2 eV. When the first bin with the 
lowest neutrino momentum states is converted to particles, only the remaining 14 bins 
are summed to get the neutrino grid Fourier modes. The phases of the neutrino particles 
are made with the same set of random numbers as were the CDM particles at z = 49. 

Second, the neutrino iV-body particle mass should correspond to the mass removed 
from the grid. 

Third, when neutrino iV-body particles are created in a certain momentum range, 
the thermal velocity they recieve lies within this range as well. This is based on the fact 
that at the conversion redshift the thermal velocity is much larger than the gravitational 
flow velocity, meaning that neutrinos in a given bin can at most have originated from 
the adjacent bins. It is therefore a very good approximation to give the neutrino 
particles a thermal velocity corresponding to the bin and assuming that there is no 
correlation between the thermal and flow components. This reduced bin mixing caused 
by the thermal velocity is fortunate, since disentangling these two components would 
be impossible. 

The thermal velocities for a given particle grid position in two succeeding 
neutrino grid-to-particle conversions are drawn with the same random numbers, but 
with momenta pointing in opposite directions. This ensures an approximate local 
conservation of momentum. In j5] local conservation was not enforced due to a higher 
neutrino particle starting redshift, and since only the average TF was used as input 
power spectrum. 

The neutrino gravitational velocity is found from the difference of two position grids 
centered around (with Aa = 0.01) our neutrino grid-to-particle conversion redshift (first 
described in [5]). When we add the gravitational flow velocity, some neutrinos will have 
a total momentum which lies outside of the bin. This is a problem when part of the 
neutrino grid is converted to particles at times when the neutrino thermal velocity is 
only a few times the gravitational velocity, though effectively being a problem only for 
the low momentum neutrinos. 

4. Results 

In Fig. [2] we show density grids from the hybrid Ai simulation (see Table [I]). Some 
features can immediately be seen from this. As expected the part of the neutrino 
distribution remaining on the grid is much less clustered on small scales than the neutrino 
particles, while the opposite is true on the largest scales. Interestingly, it can also be 
seen that the q/T = 1 bin is initially only weakly clustered on all simulated scales 
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Figure 2. Density grids from the hybrid simulation A\ with J2 m v = 1-2 eV found 
with the adaptive smoothing length kernel from [T2]. Top mosaic: z = 4. Bottom 
mosaic: z = 0. Top row: CDM, neutrino particles and neutrino grid. Middle row: 
q/T = 1, 2 and 3. Bottom row: q/T = 5, 7 and 10. Bottom mosaic: The CDM, 
neutrino particle, q/T =1,2 and 3 density grids have been raised to the power of 0.25 
to enhance the dynamical contrast. The density slices have a thickness of 20h~ 1 Mpc 
and are 512/i _1 Mpc on a side. 
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£m„ [eV] 

n„ [%] 


512 3 512 3 512 3 
512 3 512 3 
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4 oo 
1.2 1.2 1.2 
2.6 2.6 2.6 


256 3 256 3 256 3 256 3 256 3 256 3 
256 3 256 3 256 3 256 3 256 3 256 3 
4 5 6 7 10 15 
4 4 4 4 4 4 
1.2 1.2 1.2 1.2 1.2 1.2 
2.6 2.6 2.6 2.6 2.6 2.6 


256 3 256 3 256 3 256 3 
256 3 256 3 256 3 256 3 
10 10 10 10 
2 4 8 16 
1.2 1.2 1.2 1.2 
2.6 2.6 2.6 2.6 


256 3 256 3 256 3 
256 3 256 3 256 3 
10 10 10 
2 4 8 
0.6 0.6 0.6 
1.3 1.3 1.3 



Table 1. Parameters for our iV-body simulations presented in this paper. ATcdm is the 
number of CDM A^-body particles and -/V„ igr id is the size of the neutrino Fourier grid. 
Qcut/T indicates how much of neutrino momentum space is converted to particles, and 
the number of neutrino iV-body particles can be found from q cu t/T ■ N Vj g T i&. When the 
average CDM gravitational flow velocity times /flow has increased above q/T = 1, the 
conversion from grid to particles is made. £ m v is the total neutrino mass, and it is 
in all cases related to the one-particle neutrino mass, m v , by — 3m„. f2„ is the 

fraction of the critical density contributed by the neutrinos today. In all simulations 
i? BOX = 512/i _:L Mpc, and the size of the particle mesh grid is equal to the number of 
CDM particles. All simulations have a starting redshift of 49. 

(i.e. at z = 4), but quickly starts tracking the CDM component, and at z = has as 
much clustering as the next momentum bin. Note that the individual neutrino density 
grids are found from the original binning of the neutrino particles at the grid-to-particle 
conversion redshift (determined by the primordial thermal velocity), and not according 
to the actual neutrino particle velocities at redshifts 4 and 0. 

4-1. Comparing the hybrid method to its building blocks 

We have compared our new hybrid method (A\) with the full non-linear simulation (A2), 
where neutrino iV-body particles are created at z = 49, and with the simulation where 
the neutrino component stays linear on a grid (A3). 

On the left hand side of Fig. [31 it can be seen that the new method gives almost 
the same total matter power spectrum as the full non-linear simulation, except at low 
redshift where the hybrid method gives slightly less power for k > 0.2/iMpc" 1 . In the 
full non-linear simulation it can be seen from [U] that the matter power spectrum has 
not converged for 512 3 neutrino particles. The discrepancy seen in Fig. [3] is therefore 
due to a complete suppression of the noise term in the hybrid simulation, which we have 
verified by running simulations with 10 ■ 256 3 (C2) and 10 • 512 3 (Ai) neutrino iV-body 
particles with the hybrid method. 

On the right hand side of Fig. [3] we show the neutrino power spectrum for the 
different methods. With the hybrid method we are now able to simulate the non- 
linear neutrino power spectrum accurately out to k ~ 1/iMpc -1 , even with Rbox — 
512 hr 1 Mpc. The extra accuracy is of course achieved by including more neutrino 
particles, but with a comparable amount of CPU time. 
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Figure 3. Left: Difference between the linear simulation ^3 (green) and the full non- 
linear A2 (red) and hybrid A\ (blue) methods at various redshifts. Right: Absolute 
neutrino power spectra at z = 0. 



4-2. Converting part of neutrino momentum space to particles 

Since only a small fraction of the neutrino mass resides in the high momentum tail it 
is desirable to find a momentum cut-off above which the neutrinos do not contribute 
significantly to structure formation at a given scale. Such a cut-off will reduce the 
simulation time, because the timestep criterion depends in part on the total velocity of 
the particles, and particles in the highest momentum bins have thermal velocities much 
higher than the gravitional flow velocities even at z = 0. 

In Fig. (jl]) we show the total matter (left) and neutrino (right) power spectra for 
simulations with a cut-off at q C ut/T = 4, 5, 6, 7, 10 and 15. For q cut /T = 5 and 10, 86% 
and 99.7% of the neutrino mass is converted to particles, respectively. The total matter 
power spectrum is almost identical for q cut /T = 10 and 15. For q cnt /T = 5 the matter 
power spectrum increases by 2%. This is due to leakage of power across the momentum 
cut-off scale. But a 1% error in the large-scale perturbation field hardly affect the 
formation of small-scale structure. For cases (such as the study of halo profiles) where 
only the small-scale structure of neutrinos is important, q C ut/T ~ 5 would be sufficient. 

In the neutrino power spectra the same trends can be identified. There is no 
significant double counting of perturbations at small scales, since at these scales the 
linear grid does not contribute for q > q cnt . The fact that the q cu t/T < 6 simulations 
give less power at k > 1 /iMpc" 1 , can be attributed to both neutrino iV-body particle 
shot noise and a lack of non-linear corrections from the non-converted grid bins. We 
leave the relative importance of these two effects to people with too much CPU time. 

Jf.,3. The optimal grid-to-particle conversion redshift 

The optimal grid-to-particle conversion redshift is found as a trade-off between ensuring 
that non-linearities in the neutrino component is accurately simulated and that the 
finite number of neutrino iV-body particles do not introduce a significant noise term. 
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Figure 4. z = 0. Left: Hybrid simulations Bi with q cu t/T ranging from 4 to 15 
compared to q cu _t/T =15 (i^). Right: Corresponding neutrino power spectra. 
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Figure 5. *}2m u = 1.2 eV. Left: Matter power spectra from simulations Ci at z = 
with grid-to-particle conversion factors of /a ow = 2, 4, 8 and 16 relative to /fl ow = 8 
(C3). Right: Corresponding neutrino power spectra. 



The first requirement is satisfied by a high and the second by a low conversion redshift. 

In Fig. (0) we show the difference in the matter and neutrino power spectra with 
^2 m v = 1.2 eV for different conversion redshifts. The conversion criterion used was 
when the upper velocity limit of the q/T = 1 momentum bin had fallen below some 
factor /fl ow of the average CDM iV-body particle flow velocity This quantity is related 
to the strength of the gravitational field, and therefore to the ability of neutrino particles 
to cluster beyond linear order. We have run simulations with /a ow = 2, 4, 8 and 16, 
corresponding to conversion redshifts of 3.1, 5.5, 9.4 and 15.7, respectively. It can be 
seen that with /a ow = 8 the matter power spectrum can be found with 0.2% precision 
and the neutrino power spectrum with at precision at the 2-4% level. The simulation 
with /fl ow = 2 clearly lacks non-linear corrections due to the late conversion redshift. 

By running simulations with the particle grids down-scaled by 2 3 we have found 
that these results have not yet converged, though they have converged at an acceptable 
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Figure 6. J2 m v = 0.6 eV. Left: Matter power spectra from simulations Di at z = 
with grid-to-particle conversion factors of /a ow = 2, 4 and 8 relative to /a ow = 8 (-D3). 
Right: Corresponding neutrino power spectra. 



level to demonstrate the robustness of the hybrid method. Likewise, we have found that 
the lower value of the matter power spectrum with / flow = 16 compared to /a ow = 8 is 
a noise term, making the neutrino distributions decorrelate with its CDM counterpart, 
giving rise to less structure. 

In Fig. (EJ) the same case is examined with Xl m ^ = 0.6 eV, and the same general 
trends can be identified. The relative error in the matter power spectrum is proportional 
to Q u and the relative error in the neutrino power spectrum is roughly independent of 

5. Discussion and Conclusions 

We have presented a new scheme for resolving the small-scale structure of neutrinos. It 
combines features from grid simulations [S] which are fast and accurately track neutrinos 
in the linear regime with particle simulations [6] capable of resolving neutrino structures 
in the highly non-linear regime. 

The hybrid scheme is ideal for calculating the non-linear matter power spectrum for 
neutrino masses around the current upper bound '^2 i m l , ~ 0.5 eV. For lower neutrino 
masses the grid approach has an accuracy better than 1%. Furthermore, the hybrid 
scheme with separate momentum dependent initial TFs and corresponding separate 
thermal velocities reproduces our earlier results from [3] and [BJ, where the average 
neutrino TF was used. 

For Y2 m f = 1-2 eV converting 5 to 10 momentum bins at redshifts z ~ 5 — 10 
will accurately reproduce the matter power spectrum together with the neutrino power 
spectrum on small scales. For lower neutrino masses the particles can be created later 
due to their higher thermal velocities and a smaller part of the neutrino energy needs 
to be converted to particles to accurately simulate neutrino small-scale structure. 

In the present work we have focused on features in the power spectrum and limited 
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our discussion to scales of k < 1 h Mpc~ , where our code can easily keep errors under 
control at the 1% level. Achieving a comparable level of accuracy in pure particle-based 
simulations would require a prohibitive amount of CPU time. 

Even more interestingly, this hybrid scheme is capable of resolving much smaller 
neutrino structures, such as galaxy- or cluster-sized halos. This will be the topic of a 
follow-up paper [T3] . 
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